#include "2BWTSearch.h"
#include "2BWTSearch_cpu.h"
#include "mian_util.h"
#include "cuda_header.h"
#include <algorithm>
#include "coproc.h"
#include "../CompressAlign/cpp/compressTitle_cpu.h"
#include "../CompressAlign/cpp/compressOutput_cpu.h"
#include "../CompressAlign/cpp/compressOutput_gpu.h"
#include "../saIndex/alignerOutput.h"
//#include "Index.h"
//#include "PausedQueryBuffer.h"


using namespace std;

char getComplement(char base) {

        switch(base) {
                case 'a':
                        return 't';
                case 'A':
                        return 'T';
                case 't':
                        return 'a';
                case 'T':
                        return 'A';
                case 'c':
                        return 'g';
                case 'C':
                        return 'G';
                case 'g':
                        return 'c';
                case 'G':
                        return 'C';
                default: {
			printf("character : %c \n", base);
                        ERROR_EXIT("unknown base");
		}
        }
}
void reverse_complement(char* outRead, const char* inRead, const int readLength) {
        for(int i = 0; i < readLength; i++) {
                outRead[i] = getComplement(inRead[i]);
        }
        reverse(outRead, outRead + readLength);
}



void copyFileName(char * name, char * path, char * suffix)
{
	strcat(name, path);
	strcat(name,suffix);
}

char decodeInt(int i)
{
	switch(i)
	{
		case(0):
			return 'A';
			break;
		case(1):
			return'C';
			break;
		case(2):
			return'G';
			break;
		case(3):
			return'T';
			break;
	}
}

char EncodeInt(char a)
{
	switch(a)
	{
	case'A':
		return 0;
	case'C':
		return 1;
	case'G':
		return 2;
	case'T':
		return 3;
	}
	return 0;
}


uint64 reportNumMatch(const long* numMatch, const int numQuery) {
	int numDiscard = 0;
	int numAlign = 0;
	int numQueryAligned = 0;
	
	for(int i = 0; i < numQuery; i++) {
		if(numMatch[i] == -1) {
			numDiscard++;
		} else {
			numAlign += numMatch[i];
			if(numMatch[i] > 0) {
				numQueryAligned ++;
			}
		}
	}
	
	printf("reported by reportNumMatch: numDiscard: %d\n", numDiscard);
	printf("reported by reportNumMatch: numQueryAligned: %d\n", numQueryAligned);
	printf("reported by reportNumMatch: numAlign: %d\n", numAlign);

	return numAlign;
}

void LoadIndex(BWT * bwt, BWT *revbwt, char* path)
{
	int i;
	FILE *bwtCodeFile,*saValueFile, *revbwtCodeFile, *occValueFile,*revoccValueFile, *gpuOccValueFile, *gpurevOccValueFile= NULL;
	unsigned int tmp;
	unsigned int bwtCodeLengthInFile;
	const int MaxFileNameLength = 100;
//	char path[] = "/home/lumian/data/human_ref/ref_good.fa.index";
//	char path[] = "/home/lumian/data/human_ref/2bwt-index/hg19.fa.index";
	//char path[] = "data/ref.fa.index";
	char bwtCodeFileNameSuffix[] = ".bwt";
	//	char saValueFileNameSuffix[] = "/ref.fa.index.sa";
	char revbwtCodeFileNameSuffix[] = ".rev.bwt";
	char lookupTableFileNameSuffix[] = ".lkt";
	char revLookupTableFileNameSuffix[] = ".rev.lkt";
	char occValueFileNameSuffix[] = ".fmv";
	char revoccValueFileNameSuffix[] = ".rev.fmv";
	char gpuoccValueFileNameSuffix[] = ".fmv.gpu";
	char gpurevoccValueFileNameSuffix[] = ".rev.fmv.gpu";
//	char gpuoccValueFileNameSuffix[] = ".fmv.gpu";
//	char gpurevoccValueFileNameSuffix[] = ".rev.fmv.gpu";


	char bwtCodeFileName[MaxFileNameLength]="";
	char saValueFileName[MaxFileNameLength]="";
	char revbwtCodeFileName[MaxFileNameLength]="";
	char lookupTableFileName[MaxFileNameLength]="";
	char revLookupTableFileName[MaxFileNameLength]="";
	char occValueFileName[MaxFileNameLength]="";
	char revoccValueFileName[MaxFileNameLength]="";

	//char gpuOccValueFileName[] = "/home/lumian/data/human_ref/hg19_OccValueFileGPU";

	char gpuOccValueFileName[MaxFileNameLength] = "";
	char gpurevOccValueFileName[MaxFileNameLength] = "";

        //char gpuOccValueFileName[] = "OccValueFileGPU";
        //char gpurevOccValueFileName[] = "RevOccValueFileGPU";

	copyFileName(bwtCodeFileName, path, bwtCodeFileNameSuffix);
	//	copyFileName(saValueFileName, path, saValueFileNameSuffix);
	copyFileName(revbwtCodeFileName, path, revbwtCodeFileNameSuffix);
	copyFileName(lookupTableFileName, path, lookupTableFileNameSuffix);
	copyFileName(revLookupTableFileName, path, revLookupTableFileNameSuffix);
	copyFileName(occValueFileName, path, occValueFileNameSuffix);
	copyFileName(revoccValueFileName, path, revoccValueFileNameSuffix);
	copyFileName(gpuOccValueFileName, path, gpuoccValueFileNameSuffix);
	copyFileName(gpurevOccValueFileName, path, gpurevoccValueFileNameSuffix);

	unsigned int *tempBwtCode;
	unsigned int * tempPacValue;

	bwtCodeFile = (FILE*)fopen(bwtCodeFileName, "rb");
	if (bwtCodeFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open bwtCodeFile%s!\n", bwtCodeFileName);
		exit(1);
	}

	revbwtCodeFile = fopen(revbwtCodeFileName, "rb");
	if (revbwtCodeFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open revbwtCodeFile!\n");
		exit(1);
	}
#if LAYER_INDEX_2
	occValueFile = (FILE*)fopen(occValueFileName, "rb");
	if (occValueFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open occFile!\n");
		exit(1);
	}

	revoccValueFile = fopen(revoccValueFileName, "rb");
	if (revoccValueFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open revoccFile!\n");
		exit(1);
	}
#endif

#if LAYER_INDEX_1
	gpuOccValueFile = (FILE *)fopen(gpuOccValueFileName, "rb");
	gpurevOccValueFile = (FILE *)fopen(gpurevOccValueFileName,"rb");
#endif
	/*
	   if (saValueFileName != NULL && saValueFileName[0] != '\0' && saValueFileName[0] != '-') {
	   saValueFile = (FILE*)fopen(saValueFileName, "rb");
	   if (saValueFile == NULL) {
	   fprintf(stderr, "BWTLoad() : cannot open saValueFile!\n");
	   fprintf(stderr, "saValueFileName:%s\n", saValueFileName);
	   exit(1);
	   }
	   }

*/
	printf("reference file: %s\n",path);
	bwt->countValue = (unsigned int *)malloc((ALPHABET_SIZE + 1)*sizeof(int));
	fread(&bwt->inverseSa0, sizeof(unsigned int), 1, bwtCodeFile);
	fread(&bwt->countValue[1], sizeof(unsigned int), ALPHABET_SIZE, bwtCodeFile);
	Initialize(bwt, bwt->countValue[ALPHABET_SIZE]);
	bwt->countValue[0] = 0;

	bwt->bwtSizeInWord = BWTResidentSizeInWord(bwt->textLength) + WORD_BETWEEN_OCC / 2;
	bwtCodeLengthInFile = BWTFileSizeInWord(bwt->textLength);
	tempBwtCode = (unsigned int *)malloc(bwt->bwtSizeInWord * sizeof(unsigned int));
	printf("reference length: %u\n", bwt->textLength);	
	

	fread(tempBwtCode, sizeof(unsigned int), bwtCodeLengthInFile, bwtCodeFile);
	fclose(bwtCodeFile);

	DecodeBWT(bwt, tempBwtCode);
	BWTClearTrailingBwtCode(bwt);

	LoadLookupTable(bwt->lookupTable, lookupTableFileName, TABLE_SIZE);

#if LAYER_INDEX_2
	for(i=0; i<ALPHABET_SIZE+1; i++)
		fread(&tmp,sizeof(unsigned int), 1, occValueFile);
	fread(bwt->occValue, sizeof(unsigned int), bwt->occSizeInWord, occValueFile);
	fread(bwt->occMajorValue, sizeof(unsigned int), bwt->occMajorSizeInWord, occValueFile);
	fclose(occValueFile);
#endif

#if LAYER_INDEX_1
	uint32 temp;
	fread(&temp, sizeof(unsigned int), 1, gpuOccValueFile);
	for(i=0; i<ALPHABET_SIZE; i++){
		fread(&temp, sizeof(unsigned int), 1, gpuOccValueFile);
	}
	unsigned int gpuOccValueSize = (bwt->textLength + GPU_OCC_INTERVAL - 1) / GPU_OCC_INTERVAL * ALPHABET_SIZE;	
	bwt->GPUOccValue = (unsigned int *)malloc(gpuOccValueSize * sizeof(int));
	fread(bwt->GPUOccValue, sizeof(int), gpuOccValueSize, gpuOccValueFile);
	fclose(gpuOccValueFile);
	bwt->GPUOccValueInWord = gpuOccValueSize;

#endif
	
	revbwt->countValue = (unsigned int *)malloc((ALPHABET_SIZE + 1)*sizeof(int));
	fread(&revbwt->inverseSa0, sizeof(unsigned int), 1, revbwtCodeFile);
	fread(&revbwt->countValue[1], sizeof(unsigned int), ALPHABET_SIZE, revbwtCodeFile);
	Initialize(revbwt, revbwt->countValue[ALPHABET_SIZE]);
	revbwt->countValue[0] = 0;

	bwtCodeLengthInFile = BWTFileSizeInWord(revbwt->textLength);
	revbwt->bwtSizeInWord = BWTResidentSizeInWord(revbwt->textLength) + WORD_BETWEEN_OCC / 2;
	tempBwtCode = (unsigned int *)malloc(bwt->bwtSizeInWord * sizeof(unsigned int));
	fread(tempBwtCode, sizeof(unsigned int), bwtCodeLengthInFile, revbwtCodeFile);
	fclose(revbwtCodeFile);

	DecodeBWT(revbwt, tempBwtCode);
	BWTClearTrailingBwtCode(revbwt);
	LoadLookupTable(revbwt->lookupTable,revLookupTableFileName, TABLE_SIZE);

#if  LAYER_INDEX_1	
	fread(&temp, sizeof(unsigned int), 1, gpurevOccValueFile);
	for(i=0; i<ALPHABET_SIZE; i++)
		fread(&temp, sizeof(unsigned int), 1, gpurevOccValueFile);
	revbwt->GPUOccValue = (unsigned int *)malloc(gpuOccValueSize * sizeof(int));
	fread(revbwt->GPUOccValue, sizeof(int), gpuOccValueSize, gpurevOccValueFile);
	fclose(gpurevOccValueFile);
	revbwt->GPUOccValueInWord = gpuOccValueSize;
#endif
	
#if LAYER_INDEX_2
	for(i=0; i<ALPHABET_SIZE+1; i++)
		fread(&tmp,sizeof(unsigned int), 1,revoccValueFile);
	fread(revbwt->occValue, sizeof(unsigned int), revbwt->occSizeInWord, revoccValueFile);
	fread(revbwt->occMajorValue, sizeof(unsigned int), revbwt->occMajorSizeInWord, revoccValueFile);
	fclose(revoccValueFile);
#endif
}

void compressRead2Buffer(char* read, unsigned int* Queries, unsigned int readId, int queryLen)
{
	unsigned int i = readId;
	int j,k,pos = 0;
	int querySizeInWord = QUERY_SIZE_IN_WORD;

        for(j=0; j<querySizeInWord; j++) {
		for(k=0; k<CHAR_PER_WORD && pos<queryLen; k++) {
			Queries[i * querySizeInWord + j] <<= 2;
                        Queries[i * querySizeInWord + j] |= (EncodeInt(read[pos]) & 0x3);
                        pos++;
		}
                if(pos == queryLen) {
			for(; k<CHAR_PER_WORD; k++)
				Queries[i * querySizeInWord + j] <<= 2;
                }
	}
}

void addName2Buffer(char* readName, char* readNameBuffer, unsigned int* offset)
{
	int i = 0;
	while(readName[i] != '\0' && readName[i] != '\n'){
		readNameBuffer[(*offset)++] = readName[i++];
	}
	readNameBuffer[(*offset)++] = '\0';
}


bool isFASTQ(char * inputFileName) {
	bool FASTQ = false;

	if(strstr(inputFileName, "fastq") != NULL){
		FASTQ = true;
		printf("Input format: FASTQ.\n");
	} else {
		printf("Input format: FASTA.\n");
	}
	
	return FASTQ;
}

void parseInputFile(const int QUERY_NUM, FILE* fin, bool FASTQ, int* ReadLength, char* ReadNameBuffer, 
					unsigned int* ReadNum, unsigned int* Queries, 
					const bool reverse = false) {
	char c;
	unsigned int i=0, j=0, k=0, pos = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = 200;
	char ReadName[maxReadNameLength];
	char Read[maxReadLength];
	unsigned int ReadNameOffset = 0;
	char temp[200];


	fgets(ReadName, maxReadNameLength, fin);
	addName2Buffer(ReadName,ReadNameBuffer, &ReadNameOffset);
	fgets(Read, maxReadLength, fin);
	
	for(i=0; i<maxReadLength; i++) {
		if(Read[i] == '\n' || Read[i] == '\0')
			break;
	}
	(*ReadLength) = i;
	printf("%s query len: %d\n",Read, *ReadLength);
	
	int readLen = *ReadLength;
	for(j=0; j<(readLen); j++){
		c = Read[j];
		if(c == 'a' || c =='c' || c=='g' || c=='t'){
			Read[j] = toupper(c);
		}
		else if( c == 'A' || c =='C' || c=='G' || c=='T')
			Read[j] = c;
		else if(c == 'N' || c=='n')
			Read[j] = 'G';
	}
	unsigned int querySizeInWord = QUERY_SIZE_IN_WORD;
	i = 0;
	compressRead2Buffer(Read, Queries, i, readLen);
	i++;
	char* revQuery = (char*)malloc(sizeof(char)*(readLen));
	
	while(i < QUERY_NUM && !feof(fin))
	{
		fgets(ReadName,maxReadNameLength,fin);
		addName2Buffer(ReadName, ReadNameBuffer, &ReadNameOffset);
	//	printf("ReadNameBuffer: %s\n", ReadNameBuffer);
	//	printf("ReadNameOffset: %u\n", ReadNameOffset);
		
		fgets(Read, maxReadLength, fin);
	
		if(strlen(Read) != readLen + 1) {
			printf("%d different length ignore\n", strlen(Read));
			continue;
		}
		if(reverse == true) {
			reverse_complement(revQuery, Read, readLen);
			memcpy(Read, revQuery, sizeof(char)*readLen);
		}
		for(j=0; j<(readLen); j++){
			c = Read[j];
			if(c == 'a' || c =='c' || c=='g' || c=='t'){
				Read[j] = toupper(c);
			}
			else if( c == 'A' || c =='C' || c=='G' || c=='T')
				Read[j] = c;
			else if(c == 'N' || c=='n')
				Read[j] = 'G';
		}	
		
		compressRead2Buffer(Read,Queries, i, readLen);

		//pack Reads[i].query to Queries 
		if(FASTQ) {
			fgets(temp,200, fin);
			fgets(temp,200, fin);
		}
		i++;
	}
	(*ReadNum) = i;
	free(revQuery);				
}

void parseInputFile(const int QUERY_NUM, char * inputFileName, int* ReadLength, char* ReadNameBuffer, 
					unsigned int* ReadNum, unsigned int* Queries, 
					const bool reverse = false)
{
	char c;
	unsigned int i=0, j=0, k=0, pos = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = 200;
	char ReadName[maxReadNameLength];
	char Read[maxReadLength];
	FILE * fin = fopen(inputFileName, "r");
	unsigned int ReadNameOffset = 0;
	if(fin == NULL) {
		printf("can't open input file: %d\n", inputFileName);
	}
	char temp[200];
	bool FASTQ = false;

	if(strstr(inputFileName, "fastq") != NULL){
		FASTQ = true;
		printf("fastq\n");
	}

	fgets(ReadName, maxReadNameLength, fin);
	addName2Buffer(ReadName,ReadNameBuffer, &ReadNameOffset);
	fgets(Read, maxReadLength, fin);
	
	for(i=0; i<maxReadLength; i++) {
		if(Read[i] == '\n' || Read[i] == '\0')
			break;
	}
	(*ReadLength) = i;
	printf("%s query len: %d\n",Read, *ReadLength);
	
	int readLen = *ReadLength;
	for(j=0; j<(readLen); j++){
		c = Read[j];
		if(c == 'a' || c =='c' || c=='g' || c=='t'){
			Read[j] = toupper(c);
		}
		else if( c == 'A' || c =='C' || c=='G' || c=='T')
			Read[j] = c;
		else if(c == 'N' || c=='n')
			Read[j] = 'G';
	}
	unsigned int querySizeInWord = QUERY_SIZE_IN_WORD;
	i = 0;
	compressRead2Buffer(Read, Queries, i, readLen);
	i++;
	char* revQuery = (char*)malloc(sizeof(char)*(readLen));
	
	while(i < QUERY_NUM && !feof(fin))
	{
		fgets(ReadName,maxReadNameLength,fin);
		addName2Buffer(ReadName, ReadNameBuffer, &ReadNameOffset);
	//	printf("ReadNameBuffer: %s\n", ReadNameBuffer);
	//	printf("ReadNameOffset: %u\n", ReadNameOffset);
		
		fgets(Read, maxReadLength, fin);
	
		if(strlen(Read) != readLen + 1) {
			printf("%d different length ignore\n", strlen(Read));
			continue;
		}
		if(reverse == true) {
			reverse_complement(revQuery, Read, readLen);
			memcpy(Read, revQuery, sizeof(char)*readLen);
		}
		for(j=0; j<(readLen); j++){
			c = Read[j];
			if(c == 'a' || c =='c' || c=='g' || c=='t'){
				Read[j] = toupper(c);
			}
			else if( c == 'A' || c =='C' || c=='G' || c=='T')
				Read[j] = c;
			else if(c == 'N' || c=='n')
				Read[j] = 'G';
		}	
		
		compressRead2Buffer(Read,Queries, i, readLen);

		//pack Reads[i].query to Queries 
		if(FASTQ) {
			fgets(temp,200, fin);
			fgets(temp,200, fin);
		}
		i++;
	}
	(*ReadNum) = i;
	fclose(fin);	
	free(revQuery);
}

uint64 getNumAlign(uint3* sa, const uint32 numSA) {
	uint64 numAlign = 0;
	for(uint64 i = 0; i < numSA; i++) {
		numAlign += (sa[i].z - sa[i].y + 1);
	}
	
	return numAlign;
}

void sortQuery(uint3* cpuSA, unsigned int* Queries, const int numQuery, 
				BWT* bwt, BWT* revbwt) {
	preprocess(cpuSA, Queries, numQuery, bwt, revbwt);
	
	FILE* file = NULL;
	safe_fopen(&file, "query.bin", "wb");
	safe_fwrite(Queries, sizeof(unsigned int), QUERY_SIZE_IN_WORD*numQuery, file);
	safe_fclose(file);
}


void parseBinInputFile(const char* fileName, unsigned int* Queries, const int numQuery) {
	FILE* file = NULL;
	safe_fopen(&file, fileName, "rb");
	safe_fread(Queries, sizeof(unsigned int), QUERY_SIZE_IN_WORD*numQuery, file);
	safe_fclose(file);
}


void calculateNumMatch(long* numMatch, const uint3* sa, const uint32 numSA) {
	for(int i = 0; i < numSA; i++) {
		numMatch[sa[i].x] += (sa[i].z - sa[i].y + 1);
	}
}

uint32 getRefPos(uint3* SA_result, uint32 numSA, uint32* saIndex, uint32* queryId,uint32* offset, uint32* pos)
{
	uint32 i,j, k;
	k = 0;
	printf("numSA: %u\n", numSA);
	for(i=0; i<numSA; i++){
		queryId[i] = SA_result[i].x;
		offset[i] = k;
		for(j=SA_result[i].y; j<=SA_result[i].z; j++) {
			if(j > 2858012806) {
				printf("%u\n", j);
				continue;
			}
			pos[k] = saIndex[j];
			k++;
		}
	}
	printf("num alignment: %u\n", k);
	return k;
}


struct Index {
	BWT* bwt;
	BWT* revbwt;
	uint32* saIndex;
};


struct Buffer {
	//for alignment result output
	uint64 saBufSize;
	uint3* cpuSA;

	//for getRefPos
	//these three pointers shoule be allocated using allocate after the alignment process
	uint32* queryId;
	uint32* offset;
	uint32* pos;
	//these two pointers are allocated
	uint32* c_queryId;
	uint32* c_queryId_len;

	void allocate(const int numSA, const int nresult) {
		queryId =  (uint32*)malloc(sizeof(int)*numSA);
		offset =  (uint32*)malloc(sizeof(int)*numSA);
		pos =  (uint32*)malloc(sizeof(int)*nresult);
	}

	void release() {
		free(queryId);
		free(offset);
		free(pos);
	}
};


struct Read{
	char* readNameBuffer;
	unsigned int* readBuffer;
	uint32 numQuery;

	int textLength;
	int readNum;
	int maxReadNameLength;
	int readLength; 

	void writeMeta(FILE* resultFile) {
		fwrite(&textLength, sizeof(int), 1, resultFile);
		fwrite(&readNum, sizeof(int), 1, resultFile);
		fwrite(&maxReadNameLength, sizeof(int), 1, resultFile);
		fwrite(&readLength, sizeof(int), 1, resultFile);
	}
};

void runSinglePass(FILE* resultFile, Read* read, Index* index, Buffer* buf) {
	
	/*alignment*/
	uint32 numSA = 0;
	//coprocessing
	const int numCPUThread = 16;
	printf("index->bwt: %d\n", index->bwt);
	const int nresult = coproc2BWTSearch_simple(read->readBuffer, read->numQuery, index->bwt, index->revbwt, 
				numCPUThread, buf->saBufSize, buf->cpuSA, &numSA);


	/*allocate temp buffer*/
	buf->allocate(numSA, nresult);

	/*compression and write to the disk*/
	uint32 compressedSize = 0;
	uint32 compressedTitleSize = 0;
	uint32 outputSize = 0;
	read->writeMeta(resultFile);
	compressTitle(read->readNameBuffer, read->readNum, resultFile, &compressedTitleSize);
	fwrite(read->readBuffer, sizeof(int), QUERY_SIZE_IN_WORD*read->readNum, resultFile);
	sortSA(buf->cpuSA, numSA);
	getRefPos(buf->cpuSA, numSA, index->saIndex, nresult, buf->queryId, buf->offset, buf->pos);
	compressQueryId(buf->queryId, numSA, read->numQuery, buf->c_queryId, buf->c_queryId_len, &compressedSize);
	outputIdPos(buf->c_queryId, buf->c_queryId_len, compressedSize, buf->offset, 
			numSA, buf->pos, nresult, resultFile, &outputSize);

	/*free temp buffer*/
	buf->release();
}


/*
 TODO:
	parse the input read file one time.
*/
int main(int argc, char* argv[]) {


	//load the bwt index
	char* indexFileName = argv[1];
	char* queryFileName = argv[2];
	char* resultFileName = argv[3];
	char* revResultFileName = "revOutput";
	BWT *bwt, *revbwt;
	bwt = (BWT *)malloc(sizeof(BWT));
	bwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	revbwt = (BWT *)malloc(sizeof(BWT));
	revbwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	LoadIndex(bwt, revbwt, indexFileName);
	printf("Load index done. \n");
	
	uint32 numSA1 = 0, compressedTitleSize = 0;
	uint64 numAlign1, numAlign2;
	uint64 totalOutputSize = 0;
	int readLength;
	uint32 textLength = bwt->textLength;



	//load the SA index
	uint32* saIndex = NULL;
	char saFileName[200] = "";
	const char* saFileSuffix = ".sa";
	strcat(saFileName, indexFileName);
	strcat(saFileName, saFileSuffix);
	FILE* saFile = NULL;
	safe_fopen(&saFile, saFileName, "rb");
	uint32 saInterval;
	for(int i=0; i<6; i++) {
			fread(&saInterval, sizeof(int), 1, saFile);
	}
	uint32 saIndexSize = (textLength + saInterval - 1) / saInterval;
	saIndex = (uint32*)malloc(saIndexSize * sizeof(int));
	safe_fread(saIndex, sizeof(int), saIndexSize, saFile);
	safe_fclose(saFile);
	printf("Load SA done.\n");


	//init memory
	const int QUERY_NUM = 1000000;


	int i,j;
	int matchCounter = 0;
	unsigned int matchRead = 0;
	unsigned int ReadNum = 0;
	unsigned int timer = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = 100;
	const int maxReadLengthInWord = maxReadLength / 16;
	char* readNameBuffer = (char*)malloc(QUERY_NUM * maxReadNameLength);
	uint32* readBuffer = (uint32*)malloc(QUERY_NUM * 2 * sizeof(int) * maxReadLengthInWord);
	char* revReadNameBuffer = (char*)malloc(QUERY_NUM * maxReadNameLength);
	uint32* revReadBuffer = (uint32*)malloc(QUERY_NUM * sizeof(int) * maxReadLengthInWord);



	//parse the input file: TODO: parse only once!!!	
	FILE * fin = NULL;
	safe_fopen(&fin, queryFileName, "r");
	int querySizeInWord = QUERY_SIZE_IN_WORD;	
	bool fastq = isFASTQ(queryFileName);
	parseInputFile(QUERY_NUM, fin, fastq, &readLength, readNameBuffer, &ReadNum, readBuffer, false);
	rewind(fin);	
	uint32 revReadOffset = ReadNum * querySizeInWord;
	parseInputFile(QUERY_NUM, fin, fastq, &readLength, readNameBuffer, &ReadNum, readBuffer + revReadOffset, true);



	//initialize the cpu result buffer
	uint64 saBufSize = 8 *((uint64)(1024*1024*1024));
	uint3* cpuSA = (uint3*)malloc(saBufSize);


	FILE* resultFile = NULL;
	safe_fopen(&resultFile, resultFileName, "wb");
	

	const uint32 numQuery = ReadNum*2;
	Read* read = (Read*)malloc(sizeof(Read));
	Index* index = (Index*)malloc(sizeof(Index));
	Buffer* buffer = (Buffer*)malloc(sizeof(Buffer));
	read->readNameBuffer = readNameBuffer;
	read->readBuffer = readBuffer;
	read->numQuery = numQuery;
	read->textLength = textLength;
	read->readNum = ReadNum;
	read->maxReadNameLength = maxReadNameLength;
	read->readLength = readLength;
	index->bwt = bwt;
	index->revbwt = revbwt;
	index->saIndex = saIndex;
	buffer->saBufSize = saBufSize;
	buffer->cpuSA = cpuSA;
	buffer->c_queryId = (uint32*)malloc(sizeof(uint32)*numQuery);
	buffer->c_queryId_len = (uint32*)malloc(sizeof(uint32)*numQuery);


	runSinglePass(resultFile, read, index, buffer);
	
	
}

int main2(int argc, char* argv[])
{

	const int QUERY_NUM = 1000000;


	int i,j;
	int matchCounter = 0;
	unsigned int matchRead = 0;
	unsigned int ReadNum = 0;
	unsigned int timer = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = 100;
	const int maxReadLengthInWord = maxReadLength / 16;
	char* readNameBuffer = (char*)malloc(QUERY_NUM * maxReadNameLength);
	uint32* readBuffer = (uint32*)malloc(QUERY_NUM * 2 * sizeof(int) * maxReadLengthInWord);
	char* revReadNameBuffer = (char*)malloc(QUERY_NUM * maxReadNameLength);
	uint32* revReadBuffer = (uint32*)malloc(QUERY_NUM * sizeof(int) * maxReadLengthInWord);

	char* indexFileName = argv[1];
	char* queryFileName = argv[2];
	char* resultFileName = argv[3];
	char* revResultFileName = "revOutput";
	BWT *bwt, *revbwt;
	bwt = (BWT *)malloc(sizeof(BWT));
	bwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	revbwt = (BWT *)malloc(sizeof(BWT));
	revbwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	LoadIndex(bwt, revbwt, indexFileName);
	printf("Load index done. \n");
	
	uint32 numSA1 = 0, compressedTitleSize = 0;
	uint64 numAlign1, numAlign2;
	uint64 totalOutputSize = 0;
	int readLength;
	uint32 textLength = bwt->textLength;

	
	int querySizeInWord = QUERY_SIZE_IN_WORD;	
	parseInputFile(QUERY_NUM, queryFileName, &readLength, readNameBuffer, &ReadNum, readBuffer, false);
	
	uint32 revReadOffset = ReadNum * querySizeInWord;
	parseInputFile(QUERY_NUM, queryFileName, &readLength, readNameBuffer, &ReadNum, 
			readBuffer+revReadOffset, true);

	printf("query file: %s\n", queryFileName);
	printf("input #queries: %d\n", ReadNum);
	printf("query length: %d\n", readLength);
	printf("query size in word: %d\n", querySizeInWord);
	printf("readNum: %u\n", ReadNum);
	printf("MAX_NUM_ALIGN_PER_QUERY: %d\n", MAX_NUM_ALIGN_PER_QUERY);

	FILE* resultFile = fopen(resultFileName, "wb");
	if(resultFile == NULL)
		printf("cant open %s\n", resultFile);
	
	fwrite(&textLength, sizeof(int), 1, resultFile);
	fwrite(&ReadNum, sizeof(int), 1, resultFile);
	fwrite(&maxReadNameLength, sizeof(int), 1, resultFile);
	fwrite(&readLength, sizeof(int), 1, resultFile);
	
	timer = 0;
	startTimer(&timer);
	compressTitle(readNameBuffer, ReadNum, resultFile, &compressedTitleSize);
	endTimer(&timer, "compress title: ");
	totalOutputSize += compressedTitleSize;
	

	fwrite(readBuffer, sizeof(int), querySizeInWord*ReadNum, resultFile);
	totalOutputSize += querySizeInWord * ReadNum;

	uint64 nresult = 0;

//	delete(Reads);
	/////////////////////////////////////////////////////////////////////////////
	//GPU-CPU coprocessing				

        uint64 cpuBufSize = 8 *((uint64)(1024*1024*1024));
        uint3* cpuSA = (uint3*)malloc(cpuBufSize);
	const int numCPUThread = 16;	
	//=============> find all valid matches
	//coproc2BWTSearch(Queries, QUERY_NUM*2, bwt, revbwt, numCPUThread);	
	nresult += coproc2BWTSearch_simple(readBuffer, ReadNum*2, bwt, revbwt, numCPUThread, cpuBufSize, cpuSA, &numSA1);	
	//coproc2BWTSearch_simple(Queries, QUERY_NUM, bwt, revbwt, numCPUThread);
	//coproc2BWTSearch_simple(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, bwt, revbwt, numCPUThread);

	//===================> with an threshold
/*	coproc2BWTSearchLimit_simple(Queries, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM,
                                     bwt, revbwt, numCPUThread);	
*/	




	////////////////////////////////////////////////////////////////////////////
	//without CPU-GPU coprocessing
/*
	////////////////////////
	//initialization
	//CPU thread initialization

	CPU2BWTInit(bwt, revbwt);
	uint64 cpuBufSize = 4 *((uint64)(1024*1024*1024));
	uint3* cpuSA = (uint3*)malloc(cpuBufSize);
	
	//GPU thread initialization
	BWT * d_bwt = NULL;
	BWT * d_revbwt = NULL;	
	st_MemoryBuf* GPUMemBuf = NULL;
	const uint64 gpuBufSize = (2.8)*((uint64)(1024*1024*1024));
	GPU2BWTInit(&d_bwt, &d_revbwt, bwt, revbwt, &GPUMemBuf, gpuBufSize);
				

	//processing either by GPU or CPU
	//GPU===>
	timer = 0;
	startTimer(&timer);
//	#ifdef GPU_NO_THRESHOLD 
	//======> find all valid alignments
	nresult += GPU2BWTSearch_shareBuf(readBuffer, ReadNum*2, d_bwt, d_revbwt, GPUMemBuf, cpuSA, &numSA1);
//	nresult += GPU2BWTSearch(Queries, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//	nresult += GPU2BWTSearch_noConcurent(Queries, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//	nresult += GPU2BWTSearch1Kernel(Queries, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//       GPUMemBuf->reset();
//  	nresult += GPU2BWTSearch_shareBuf(revReadBuffer, ReadNum, d_bwt, d_revbwt, GPUMemBuf, cpuSA, &numSA1);
//  	nresult += GPU2BWTSearch(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//   	nresult += GPU2BWTSearch_noConcurent(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
//   	nresult += GPU2BWTSearch1Kernel(Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
        GPUMemBuf->reset();
	cudaThreadSynchronize();
	printf("\n\n===============> #alignment: %lu\n", nresult);
	endTimer(&timer, "GPU processing");
*/




	//=======> with a threshold
/*	long* d_numAlign = NULL;
	GPUMALLOC((void**)&d_numAlign, sizeof(long)*QUERY_NUM);	
	cutilSafeCall(cudaMemset(d_numAlign, 0, sizeof(long)*QUERY_NUM));
        GPU2BWTSearchTwoStrandLimit(d_numAlign, Queries, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, d_bwt, d_revbwt, GPUMemBuf);
        cudaThreadSynchronize();
	GPUMemBuf->close();
	long* gpu_numAlign = (long*)malloc(sizeof(long)*QUERY_NUM);
	FROMGPU(gpu_numAlign, d_numAlign, sizeof(long)*QUERY_NUM);
	reportNumMatch(gpu_numAlign, QUERY_NUM);
*/
/*
	//CPU===>
	long* numMatch = (long*)malloc(sizeof(long)* ReadNum);
	memset(numMatch, 0, sizeof(long)* ReadNum);
	//======>find all valid
	timer = 0;
	startTimer(&timer);
//	numSA1 = CPU2BWTSearch(cpuSA, readBuffer, ReadNum*2, bwt, revbwt);
	
//	fwrite(&numSA1, sizeof(int), 1, resultFile);
			


//	CPUPrintResult(cpuSA, numSA1, resultFileName);
//	numAlign1 = reportCPUResult(cpuSA, numSA1);
//	numSA2 = CPU2BWTSearch(cpuSA, revReadBuffer, ReadNum, readLength, bwt, revbwt);
//	CPUPrintResult(cpuSA, numSA2, revResultFileName);

//	numAlign2 = reportCPUResult(cpuSA, numSA2);
	printf("\n\n===============> #alignment (all valid hits): %lu\n", numSA1);
//	cutilSafeCall(cudaThreadSynchronize());
	endTimer(&timer, "CPU all matches");
*/
	printf("\n\n\n");

/*	//======>with the #match threshold
	printTitle("CPU2BWTSearchTwoStrandLimit");
	timer = 0;
	startTimer(&timer);
	CPU2BWTSearchTwoStrandLimit(numMatch, cpuSA, Queries, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM,
                                        QUERY_NUM, bwt, revbwt);
	cutilSafeCall(cudaThreadSynchronize());
	endTimer(&timer, "CPU 2BWTSearchLimit");
	reportNumMatch(numMatch, QUERY_NUM);
*/

/*
	//===>check result for the limit implmentation
	long* numMatch_gold = (long*)malloc(sizeof(long)*QUERY_NUM);
	memset(numMatch_gold, 0, sizeof(long)*QUERY_NUM);
        numSA1 = CPU2BWTSearch(cpuSA, Queries, QUERY_NUM, bwt, revbwt);
	calculateNumMatch(numMatch_gold, cpuSA, numSA1);	
        numSA2 = CPU2BWTSearch(cpuSA, Queries + QUERY_SIZE_IN_WORD*QUERY_NUM, QUERY_NUM, bwt, revbwt);
	calculateNumMatch(numMatch_gold, cpuSA, numSA2);
	uint64 gold_numAlign = 0;
	for(int i = 0; i < QUERY_NUM; i++) {
		if(numMatch_gold[i] <= MAX_NUM_ALIGN_PER_QUERY) {
			gold_numAlign += numMatch_gold[i];
		}
	}
	printf("numMatch_gold with the threshold %d: %lu\n", MAX_NUM_ALIGN_PER_QUERY, gold_numAlign);

*/

	//////////////////////////////////////////////////////////////////////////////
	 
	printf("numSA: %u\n", numSA1);
	sortSA(cpuSA, numSA1);
	printf("sorting done\n");

#ifdef __DEBUG__
	uint32 numAlign = 0;
	for(i=0; i<numSA1; i++) {
		numAlign += cpuSA[i].z - cpuSA[i].y + 1;
	}
	nresult = numAlign;
	printf("numAlign: %u\n", numAlign);


	FILE* saResultFile = fopen("saInterval", "wb");
	if(saResultFile == NULL) {
		printf("can't open saInterval file\n");
		return 0;
	}
	fwrite(&numSA1, sizeof(int), 1, saResultFile);
	fwrite(&nresult, sizeof(int), 1, saResultFile);
	fwrite(&textLength, sizeof(int), 1, saResultFile);
	fwrite(cpuSA, sizeof(uint3), numSA1, saResultFile);
	fclose(saResultFile);

	printf("free\n");
	//GPUBWT->close();
#endif
	startTimer(&timer);	

	free(bwt->lookupTable->table);
	free(revbwt->lookupTable->table);
	free(bwt->GPUOccValue);
	free(revbwt->GPUOccValue);
	free(bwt);
	free(revbwt);
	
	free(readNameBuffer);
	free(readBuffer);
	free(revReadNameBuffer);
	free(revReadBuffer);

	endTimer(&timer, "free bwt index ");

	uint32 *queryId, *offset, *pos, *c_queryId, *c_queryId_len, compressedSize;
	uint32 numQuery = ReadNum * 2;


	c_queryId = (uint32*)malloc(numQuery * sizeof(int));
	c_queryId_len = (uint32*)malloc(numQuery * sizeof(int));
	
	printf("\n\n\n");
	startTimer(&timer);	
	getRefPosResult(cpuSA, numSA1, nresult, textLength, &queryId, &offset, &pos, indexFileName);
	endTimer(&timer, "time for get ref pos:  ");
	
	printf("\n\n\n");
	startTimer(&timer);	
	compressQueryId(queryId, numSA1, numQuery, c_queryId, c_queryId_len, &compressedSize);
	endTimer(&timer, "time for compress queryId: ");
	
	printf("\n\n\n");
	startTimer(&timer);	
	uint32 outputSize = 0;
	outputIdPos(c_queryId, c_queryId_len, compressedSize, offset, numSA1, pos, nresult, resultFile, &outputSize);
	totalOutputSize += outputSize;
	endTimer(&timer, "time for write output: ");
	
	printf("\n\n\n");
	printf("total compressed output size: %u\n",  totalOutputSize);
	
	free(c_queryId);
	free(c_queryId_len);

/*	
	FILE* saResultFin = fopen("saInterval", "rb");
	fread(cpuSA, sizeof(uint3), numSA1, saResultFin);
	fclose(saResultFile);

	printf("numSA: %u\n", numSA1);
	char* saFileName = "/home/luo_group/data/human_ref/ref_good/ref_good.fa.index.sa";
	FILE* saFile = fopen(saFileName, "rb");
	uint32 saInterval;
	for(i=0; i<6; i++)
		fread(&saInterval, sizeof(int), 1, saFile);

	uint32 saIndexSize = (textLength + saInterval - 1) / saInterval;	
	uint32* saIndex = (uint32*)malloc(saIndexSize * sizeof(int));
	uint32* queryId =  (uint32*)malloc(sizeof(int) * nresult);
	uint32* offset =  (uint32*)malloc(sizeof(int) * numSA1);
	uint32* pos =  (uint32*)malloc(sizeof(int) * numSA1);

	startTimer(&timer);
	fread(saIndex, sizeof(int), saIndexSize, saFile);
	endTimer(&timer, "read sa");
	fclose(saFile);

	startTimer(&timer);
	getRefPos(cpuSA, numSA1, saIndex, queryId, offset, pos);
	endTimer(&timer, "get ref pos ");
	free(saIndex);
	
	startTimer(&timer);
	fwrite(&nresult, sizeof(int), 1, resultFile);
	fwrite(&numSA1, sizeof(int), 1, resultFile);
	fwrite(queryId, sizeof(int), numSA1, resultFile);
	fwrite(offset, sizeof(int), numSA1, resultFile);
	fwrite(pos, sizeof(int), nresult, resultFile);
	endTimer(&timer, "write queryId and pos to disk ");
		
*/
/*	
	FILE * tempFout = fopen("out_SAResult", "w");
	for(i=0; i<numSA1; i++){
		fprintf(tempFout,"%u %u %d %d %d %d %d\n", cpuSA[i].queryId, cpuSA[i].pos, cpuSA[i].variation[0], 
			cpuSA[i].variation[1], cpuSA[i].variation[2], cpuSA[i].variation[3], cpuSA[i].variation[4]);
	}
	fclose(tempFout);

*/

/*
	Info * info = (Info*)malloc(sizeof(info));
	info.numMismatch  = (uint8* ) malloc(numSA1);
	info.misPos[0] = (uint8* )malloc(numSA1);
	info.misPos[1] = (uint8* )malloc(numSA1);
	info.misBase[0] = (uint8* )malloc(numSA1);
	info.misBase[1] = (uint8* )malloc(numSA1);

	info.readStringLength = readLength * readNum;
	info.refStringLength = textLength;

	findMisMatch(readBuffer, ref, readLength, queryId, pos, numSA1, info);


	uint64 curPos = ftell(resultFile);

//	test_compress_blockSize(cpuSA, numSA1, resultFile, 1000000, curPos);
//	test_compress_blockSize(cpuSA, numSA1,resultFile, 10000000, curPos);
//	test_compress_blockSize(cpuSA, numSA1,resultFile, 100000000, curPos);	
	printf("here\n");	
//	unsigned int timer = 0;
	startTimer(&timer);
	test_compress_blockSize(queryId, pos, info, numSA1, resultFile, numSA1, curPos);	
	endTimer(&timer, "total cpu ");
	printf("\n\n\n");
	
	startTimer(&timer);
	gpu_test_compress_blockSize(queryId, pos, info, cpuSA, numSA1, resultFile, numSA1, curPos);	
	endTimer(&timer, "total gpu");

	free(info.misPos[0]);
	free(info.misPos[1]);
	free(info.misBase[0]);
	free(info.misBase[1]);
	free(info.numMismatch);
*/
	free(cpuSA);
/*	free(pos);
	free(offset);
	free(queryId);
*/	fclose(resultFile);
	return 0;
}
